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Abstract Observations of the interplanetary shock provide us with strong evidence of 
particle acceleration to multi-MeV energies, even up to GeV energy, in a solar flare or 
coronal mass ejection (CME). Diffusive shock acceleration is an efficient mechanism for 
particle acceleration. For investigating the shock structure, the energy injection and en- 
ergy spectrum of a CME-driven shock, we perform dynamical Monte Carlo simulation of 
the 14-Dec-2006 CME-driven shock using an anisotropic scattering law. The simulated 
results of the shock fine structure, particle injection, and energy spectrum are presented. 
We find that our simulation results give a good fit to the observations from multiple space- 
craft. 

Key words: Acceleration of Particles, CME-driven Shock, Solar Energetic Particles, 
Numerical Simulation 



1 INTRODUCTION 

It is widely accepted that there are two classes of solar energetic particle (SEP) events, although recen t 
observations indicate that the actual processes may be much more complicated dPick & Vilmerl 20081) . 



The f irst class is normal impulsive SEP events, which are connected with the large solar flare dMillerl 
119971) . The second class is gradual SEP events, which are re sponsible by diffusive shock acceleration 
(PSA) as s ociati ng with fast coronal mass ejections (CMEs) dCane. Reames. & von Rosenvingel 1 1 99 lb 
lYan et alll2006h . In solar magnetic connecti on region, CME and flares are two type of m anifestations 
of the same magnetic energy release process dWang et aUll996tlZhang et al.Ll200lil2007l) . Both CMEs 
and flares result in particle acceleratio n that constitute an SEP event. But which manifestation dominates 
the particle injection is still not clear dLi et al.l 120091; iLe et all 1201 \i loin & Shalchil l2009t) . Some nu- 
merical models suggest that mixed particle acceleration by both flares and CME-driven shocks provide 
much better fits to the in-situ observations. Since the particle injection process is connected with the 
complicated nonlinear effects in the particle acceleration processes and also there exists difference of 
injection mechanism in two type of manifestations, so we just put forward to a pure shock numerical 
model to calculate the particle injection in the CME-driven shock. We expect that it would be helpful 

for understanding the particle injection problem in the SEP events. 

PSA th e ory was first introduc e d in t he later of 1970's ( iKrvmskvl dl977l) : lAxford et all dl977l) ; 
iBelll d 1978b ; iBlandford & Ostrikerl dl978l) . In the past several decades, the accumulation of the in- 
creasingly observational data from many spacecraft investigated the nonlinear diffusive shock ac- 
celeration (NLPSA) mechan ism, which is the most efficient accelerator in many astrophysical and 
space physical environments dBednarz & Ostrowskil 19991; Malkov & Prurvl 1200 it iBvkov et al.l 120091 ; 



iBvkov & Treumannl 120111 iLu. Xia , & Wangl 120061: IZhang. Bi. & Hul |2006[) . With the development of 
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technology in observational equipments, especially in spacecraft working in deep space, there are a lot 
of models for modeling the various nonlinear interaction with the diffusive shoc k acceleration. Severa l 
main approaches for studying the n onlinear DSA includes: t he two-fluid mode l (iDrurv & Volki 1981 



Dmrv, Axford, 



2007; Zirakashvili, 



& Sum mers 
2007t 



19821); the numeric a l mod el dBerezhko & Volki l2000t iKang & Jones , 



Verkhoglvadova et all 20101): the s t ationary or dynamical Mont e Carl o 
model dEllison & Eichlerl 1 1984 iKnerr. Jokipii & Ellisonl 1 19961; IVladimiroy. Ellison & Bvkovi 120061) : 
si dj 



20001; ICaprioli. Amato & Blasil |2010|) and etc. Among 
Carlo method addresses the nonlinear effects of DSA by assum- 
a random walk under a certain scattering law 



the semi-analytical model dMalkov et alj 
these approaches, the Monte 
ing that the entire particle population undergoes 

dEllison. Mobius & Paschmannlll990l:lKnerr. Jokipii & Ellisonl 1 1996tlWang & Yan[l2oTlt) . 

There are three important non-linear processes of DSA t heory including t he particle injection, par- 
ticle confinement, and shock robustness dMalkov & Drurvl 1200 lb iHuL l2009h . Owing to the fact that 
walking processes of the particles can be controlled self-consistently in the Monte Carlo method, the 
Monte Carlo method has an advantage for simulating the particle injection. We have already studied the 
energy translation processes of an Earth-bow shock using the d ynamical Monte Ca rlo method with a 
prescribed multiple anisotropic scattering angular distributions ( Wang & Yard, 120121) . We find that the 
acceleration efficiency increases as the dispersion of the scattering angular distribution increases from 
an anisotropic case to an isotropic case. Here, we will further investigate this important particle injection 
problem in the CME-driven shock using the dynamical Monte Carlo method. There exist a few differ- 
ent properties between the Earth-bow shock and the CME-driven shock: Firstly, Earth-bow shock has 
a stationary downstream bulk flow but the CME-driven shock has a dynamical downstream bulk flow; 
Secondly, the CME-driven shock front has an opposite motion compared with the Earth-bow shock's 
evolution; Thirdly, the CME-driven shock has an extended plane shock front structure near the Earth, 
but the Earth-bow shock front has a stationary bow shock geometry. We predict those differences would 
produce a different non-linear properties including the evolution of the shock fine structure, energy in- 
jection rate and even the energy spectral shape. This paper will focus on the understanding some of 
the non-linear properties of the planetary CME-driven shock. The 14-Dec-2006 shock event was fortu- 
itous as it provides us an opportuni ty for applying the d ynamical Monte Carlo package code, which was 
developed on the Matlab platform dWang & Yanll201 ll) . 

This paper is structured as follows: In Section [2] we present the specific observations for the 14- 
Dec-2006 CME-driven shock event; The detailed description of the method is given in Section [3j We 
present the simulated results and discussions in SectionlU Finally, Section[5]presents the summary and 
some conclusions. 



2 OBSERVATIONS 

The unusual group of CME-driven shock events of solar cycle 23 was observed in December 2006 at 
the solar active region 10930. Halo CMEs were observed by the LASCO coronagraphs in association 
with the events of 13 and 14 December, with speeds of 1774km/s and 1042km/s, respectively. Because 
the 14 December solar event was better magnetically connected to the Earth, so it provided the best 
opportunity for testing the nonlinear effect and efficiency of the diffusive shock acceleration (DSA) 
mechanism. As shown in FigureQ] an overview of key parameter observations from the Proton Monitor 
(PM) instruments on Wind/SWE for the CME shock event of 14 December 2006 is given in detail. This 
event originated on the western hemisphere of the Sun. It showed an abrupt fluctuation in intensity of 
proton density and solar wind thermal speed during the decay of the 1 3 December solar event. The initial 
particle increase following the 14 December solar event was seen in the higher energy range, as expected 
for velocity dispersion. There was also a hi gher background at the lower energy associated with the 13 
December solar event and the related shock dvon Rosenvinge et afll2009h . Simultaneously, many spikes 
were also detected to be superposed on the radio continuum in the frequency range 2.6-3.8 GHz by the 
digital spectrometers of NAOC. These spikes were found to have complex str uctures associated with 
other radio burst signatures connecting with the in-situ SEP event observations dWang etalll2008l) . 
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Fig. 1 The plot shows the key parameters of the 14-Dec-06 shock event in the Wind space- 
craft, and the data come from http : / / cdaweb.gsfc.nasa.gov/cdaweb. 



Both Wind and ACE were in orbit around the LI Lagrangian point ~ 1.5 million km upstream 
of the Earth. Simi l ar inte nsity modulations were observed at Wind and ACE. As has also been noted 
by iMulligan et al.l d2008l) . the variations of the particle intensity and smooth magnetic fields observed 
by near-Earth spacecraft occurred in the duration of the interplanetary coronal mass ejection (ICME) 
driving the shock on 14 December which was related to the 13 December solar event. Solar wind ob- 
servations from ACE show evidence of the presence of the ICME, which had the enhanced magnetic 
field and a smooth rotated "magnetic cloud" in the upstream shock and the intervening sheath region, 
respectively. According to the Wind magnetic clo ud list, the axis orientation of this "magnetic cloud" 
was 9 = 27°, (j) = 85°. In addition. lLiu et al.l(l2008l) estimate that the "cloud" axis direction was =—57°, 
<f> = 81° in GSE coordinates. Thus, both agree that the axis was closely aligned west to east but differ in 
whether it was inclined north or south, most likely because different intervals were considered in their 
analysis. 
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Fig. 2 The upper panel represents the proton density profile vs its position at the end of the 
simulation. The middle panel denotes the solar wind thermal velocity profile in the local frame 
vs the time. The lower panel indicates the bulk flow speed profile vs its position at the end 
of the simulation. The vertical lines in the upper and lower panels both show the final FEB 
position at the end of simulation. 



3 THE METHOD 
3.1 Physical model 

We consider a plane-parallel shock where the supersonic flow moves from the Sun to the Earth (in the 
rest frame) along the x-axis direction. The shock was observed by Wind, SOHO, and ACE spacecraft 
near the Earth in the location of the first Lagrangian point LI ~ 1.5 million km (~ 250i? e , where R e is 
the radius of the Earth) upstream of the Earth on 14 December. All trajectories of the spacecraft in the 
348th day corresponding to the 14 December 2006 are shown in Figure[3] With the CME-driven shock 
propagating from the Sun along x-axis to the Earth, its shock front were encountered by Wind, SOHO, 
and ACE spacecraft located in Xqse between the 250i? e and 180i? e upstream to the Earth. These three 
spacecraft moved about 10i? e distance in their obits on the 348th day. The distances of all these three 
spacecraft from the Sun-Earth line were within 50Re along the Yqse an d Zqse directions. The 14-Dec- 
2006 shock event originated from the western hemisphere of the Sun with an interplanetary " magnetic 
cloud" axis orientation of 9 = 27°, <f> = 85°. And the actual trajectory of Wind spacecraft at that moment 
is just tangent to the Sun-Earth line with an angle <fr = 80° as shown in the lower panel of the Figure [3] 
As far as the position of Wind spacecraft is concerned, the observed CME-driven shock is just a parallel 
diffusive shock. So the observation of Wind spacecraft provided an example of semi-parallel shock for 
applying our dynamical Monte Carlo code to understand the particle injection problem of DSA theory. 
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Table 1 The Parameters of the Simulated Shock 



Simulated Parameters 


UllllCIlMUIllCsa Veil LIC 


OUcllCQ ValLlC 


Physical 


Upstream bulk speed 
Downstream bulk speed 
Relative inflow velocity 
Inflow sonic Mach number 
Thermal speed 
Scattering time 


u u =-0.4467 
u d =-0.7742 
Au=0.3275 

M=17.5 
u =0.0342 
r=0.3333 


-600km/s 
-1042km/s 
442km/s 

46km/s 
0.052s 


Numerical 


Box size 
Total time 
Time step size 
Number of zones 
Initial particles per cell 
FEB distance 


T max =24W 

dt=mo 

nx=6W 
n =650 
X Jeb =9Q 


6.3minutes 
0.0053s 

3-R e 



Notes: The physical parameters are taken from the Wind spacecraft, and the numerical parameters are decided by 
the 14-Dec-2006 CME-driven shock. 

The important physical parameters of this simulation include the upstream bulk flow velocity (u u ), 
the downstream bulk flow velocity (ud), the relative bulk inflow velocity difference (Am), the inflow 
sonic Mach number (\ud\/c s ), which is 17.5 (where c s = {^kT/m) 1 / 2 , c s is the upstream sound speed), 
the upstream thermal velocity [vq = (kT/m) 1 / 2 ], and the constant scattering time (r), which is 2/5 
times of the scattering time (tq = 0.13s) used by Knerr, Jokipii & Ellison (1996) in the Earth-bow 
shock simulation. Since there are some differences in the shock geometry between the CME-driven 
shock (i.e. at LI point) and the Earth-bow shock, we chose the scattering time 0.4 times smaller than 
that in the Earth-bow shock, which is equivalent to a 2.5 times larger FEB distance than that in the 
Earth-bow shock. The specific physical parameters and numerical parameters are listed in the Table Q] 
with their dimensionless values and scaled values, respectively. 

3.2 Mathematical model 

According to the observed 14-Dec-2006 CME-driven shock, the schematic diagram of the simulation 
box can be designed as one-dimensional parallel shock along the x-axis direction. As shown in Figure 
HJ the initial particles with a relative bulk flow speed difference (Au) move from the right to the left. 
The initial particles have a background Maxwellian thermal distribution with an initial temperature 
(To = mv 2 /k) in the local frame. To begin and maintain the shock simulation, particles are assumed 
to flow into the simulation box from the pre-inflow box (PIB) at the right boundary. Then, with the 
continuous particle flow moving forward one time step, only those particles which move into the main 
simulation box are actually added to the simulation. This process naturally leads to a flux-weighted 
inflow population. At the left boundary of the box, a reflective wall acts to produce a CME-driven shock 
moving from the left to the right. Considering the geometry of the 14 December shock event, we just 
follow the parallel component of the CME-driven shock observed in Wind spacecraft. 

Figure|4]also shows one typical particle and its local (Vl) or box frame (V) velocity in the upstream 
region and downstream region, respectively. The majority of the incoming particles cross the shock 
front only once from the upstream region to the downstream region and stay in the downstream region. 
A small portion of the particles can effectively scatter off the resonant MHD wave self-generated by 
the energet ic particles and return to the up stream region to re-cross the shock front for additional en- 
ergy gains dLiu. Petrosian & Masonll2004 . Thus an anisotropic energetic particle distribution, but not a 
strict Maxwellian distribution, is produced in the diffusive regions. It is this elastic interaction between 
individual particles and the collective background that allows the Fermi acceleration to occur. 

The position of the FEB could coincide with a location upstream of the shock where particles are 
no longer able to scatter effectively and return to the shock. A reasonable FEB farther out in front of 
the shock moves, companying with the shock front, at the same shock evolutional velocity V s h- This 
constant FEB distance is acted to inform a precursor region which is showed by the shadow area in the 
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Fig. 3 The diagrams show the realistic obits of the near-Earth spacecraft. The obit data are 
taken from http : / / cdaweb.gsfc.nasa.gov/cdaweb. 



middle of the Figure[4] If one particle archives to the highest energy, and exceeds the position of the FEB 
in front of the shock, it will be taken as the escaped particle and removed from the simulation system. 
According to the actual motion of Wind spacecraft in the duration, the spacecraft moved about 10i? e 
distance in their obits on the 348th day. To simulate the shock formation and evolution, the total length 
of the simulation box is set to be 10i? e , the length of the FEB is set to be ~3 R P . The scattering time is 
set to be r = 0.052s on the basis of the Earth-bow shock model dKnerr. Jokipii & ElUsonl[T996h . 

The important numerical parameters include the box size (x max ), the time to evolve the whole sys- 
tem (t max ), the number of grid zones (n x ), the initial number of particles per zone (no), and the size 
of the time step (dt). Because of the similar character of the plasma flow near the Earth, we take some 
numerical parameters as in the Earth-bow shock model dKnerr. Jokipii & Ellisonl Il996l) . Specifically, 
the total box length x max = 300 is divided into n x = 600 grids, with each grid length being Ax = 1 /2; 
the total time t max = 2400 is divided into n t = 72000 steps by dt, with each step being dt = 1/30. 
All numerical parameters are listed in Table Q] The physical parameters and the numerical parameters 
constitute the whole simulation parameter list. As shown in Table Q] each dimensionless value is corre- 
sponding to its scaled value. The scale factors for distance, velocity, and time are x sca i e = 10i? e /300, 
Vscaie = 442fcms _1 /0.3275, and t scaie = x sca i e /v sca i e , respectively. 
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The Schematic Diagram of the Simulation Box 
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Fig. 4 The schematic diagram of the simulation box. The left reflective wall acts ICME to 
produce the ICME-driven shock prorogating from the left boundary to the right boundary of 
the box. 



The presented simulation s apply the same steps like the Earth-bow shock model 
dKnerr. Jokipii & EllisonL Il996h including three sub steps: (i) All the particles moving with their 
velocities in the simulation box along the x axis direction, (ii) Summing particle masses and velocities 
over each background computational grid, (iii) Invoking the scattering angular distribution law. The 
particle diffusive processes in the presented simulations are dominated by the Gaussian scattering 
angular distributions. The scattering rate is R s = dt/r, which implies that only this fraction of particles 
is able to scatter off the scattering center frozen in the background fluid. The candidate does not change 
its route until it is selected to scatter once again. So the particle's mean free path is proportional to the 
local thermal velocities in the local frame with 



A = V L ■ t. 



(1) 



For an individual proton, the grid-based scattering center can be seen as a sum of individual momenta. 
So these scattering processes can be taken as the elastic collisions. In an increment of time, once all 
of the candidates complete these elastic collisions, the momentum of the grid-based scattering center 
is changed. In turn, the momentum of the grid-based scattering center will affect the momenta of the 
individual particles in their corresponding grid in the next increment time. One complete time step 
consists of the above three substeps. The total simulation temporally evolves forward by repeating this 
time step sequence. To calculate the scattering processes accurately and produce an exponential mean 
free path distribution, the time step should be less than the scattering time (i.e. dt < r). 

The scattering angles consist of two variables: 89 and 6<f>. Once a particle has a collision with 
the background scattering centers, its pitch angle becomes 9'=9+89, and the azimuthal angle becomes 
4>'=4>+S(j>, where 56 is the variation in the pitch angle 6, and 8(f) is the variation in the azimuthal angle <f>. 
The pitch angles 9 and 9' are both in the range < 9, 9' < ir, and azimuthal angles cf> and cf>' are both in 
the range < </>, <f>' < 2ir on the unit sphere. The variation in the pitch angle 89 and azimuthal angle 8<f> 
are composed of the scattering angle, and its anisotropic character is described by the Gaussian function 
f(86, 8(f>). Here, we will just present the results of the CME-driven shock using the Gaussian scattering 
angular distribution with a standard deviation value of a = it. 
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4 RESULTS 
4.1 Data analysis 

Since the individual particle energy can be examined at any given time in the simulation, so the energy 
function over time can be obtained. At first, we can calculate the necessary energy functions for further 
analysis. In this simulation, we obtained the total energy function in the box, the loss energy function 
escaped from the FEB, and the injected energy function which is the energy summation of the injected 
energetic particles from the downstream region at the local velocity of Vl = Uq over time. Then, at the 
end of the simulation, we obtained the final values of the total energy E to t = 3.5666, the energy loss 
Ei oss = 0.2010, and the energy injection E in j = 0.5464, respectively. The final energy injection rate 
Rinj, which represents the acceleration efficiency, can be defined by the formula as follows. 

Rinj — Ei n j I E to t (2) 

The injection rate is so important for a CME-driven shock, because it is connected with the facts that the 
shock how distribute itself energy to accelerate cosmic ray (CR) and to "heat" the thermal background 
plasma. By a series of simulations, we give a plausible injection rate with a value of Ri n j = 15.32% for 
the 14-Dec-2006 shock. Under this condition, we obtain the maximum energy particle with the dimen- 
sionless value of VL max = 20.2609 and the scaled value of E max = ?>.%&%AMeV. In addition, because 
there exists some energy losses in the simulation system, the shock fine structures do not completely 
agree with the situ observations. Finally, according to the DSA theory, the energy spectrum index can 
be calculated based on the simulated compression ratio (i.e., T = (r + 2)/[2(r — 1)]). We calculated 
the total energy spectral index with a value of I\ ot = 0.8406 and the vicious subshock energy spectral 
index with a value of T su i, = 1.1074, respectively. 

Figure [5] shows the simulated energy spectra. The first plot shows the energy spectrum with the 
"double-peak" structures averaged over the entire simulation box at the end of the simulation. The 
second plot shows an energy spectrum with a "power-law" tail averaged only over the downstream 
region at the end of the simulation. In this plot, the thick solid curve with a narrow peak represents 
the initial Maxwellian thermal energy spectrum in the shock frame. As viewed from the first plot, the 
double peaks imply that there exist two thermal particle distributions in the entire simulation box: the left 
peak represents the "heated" downstream flow distribution and the right peak represents the Maxwellian 
distribution in the unshocked upstream flow. Turn to look at the second plot, we find the final extend 
energy spectrum at the left of the panel shows several decade times wider than that in the initial energy 
spectrum at the right of the panel. This means that there exist a large temperature difference between 
the shocked downstream region and the unshocked upstream region. 

As shown in Figure [6] the spectra of protons in the two largest December 2006 SEP events by 
ACE, STEROEO, and SAMPEX instruments are reported. The particle intensities started increasing 
at the beginning of the December 5 and the other at the onset of the December 13 event. These two 
events both have spec t ra tha t roll over in a similar fashion beyond ~50MeV, a s in the 20-Jan-2005 SEP 
event dMewaldt et all l2008t IWang. Zhao. & Zhoul l2009t iBartoli et all l2012t IWang et all 120101) . The 
fitting energy spectral shape of the 12-14 December 2006 events are showed in the lower panel and 
the spectral index is marked as a value of E~ im in the lower energy range. The predicted subshock 
energy spectral index (T su b = 1.1074) from our simulation is consistent with the observed energy 
spectral index in the lower energy range. Owing to computer constraints on the size of the simulation 
grid, this simulated energy spectrum is just in the range from ke V to MeV. We speculate that the second 
"roll over" on the higher energy spectrum could be obtained if a larger simulation box size is used. 
This will be investigated in a future simulation. There are two conditions suggesting that the "roll-over" 
would be reproduced at high energy: (1) The FEB distance decides the maximum diffusive length (i.e. 
FEB = X max = t ■ Pmax)- If we enlarge the FEB distance and the total simulation box, we can 
obtain the larger P max in the new simulation system. (2) In the Figure|6l we can see the first power law 
E~~ lm as the input function of the second power law E~ 2Ab at the high energy range. Simultaneously, 
in the Figure |5J we can see the heated Maxwellian thermal distribution, which would be represented by 
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Fig. 5 The extend energy spectra in two plots are calculated over the entire simulation region 
and the only downstream region at the end of the simulation, respectively. The solid extended 
curve with a "power-law" tail in each plot represents the final shocked energy spectrum. The 
thick curve with a narrow peak denotes the initial Maxwellian energy spectrum. 



a similar power law E 5 averaged over the respective energy range, as the input function of the first 
power law E" 1 - 1074 . 

4.2 Shock structures 

At the end of the simulation, the simulated shock with the specific parameter values are given as fol- 
lows: the shock position X s h = 121.5, the FEB position Xf = 31.5, the shock evolutional velocity 
V s h = —0.0744, the subshock velocity V su b = 0.2103, total compression ratio r to t = 5.4034, sub- 
shock's compression ratio r su b = 3.4697, total energy spectral index T to t = 0.8406, subshock's energy 
spectral index T su i, = 1.1074, particle injection rate R in j = 15.32%, energy loss Ei oss = 0.2010, 
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Fig. 6 Fluency spectra of the protons measured in the two largest December 2006 SEP events 
by multiple spacecraft. The energy range is from 5 to 100 MeV adapted from Mewaldt et al. 
(2008). 



the maximum energy particle's local velocity VL max = 20.2609, and the maximum energy particle 
E max = 3.8684MeV. 

First, we present the entire shock evolution with the temperature profile of the time sequences as 
shown in Figure [7] The supersonic continuous inflow with an initial Maxwellian thermal velocity vq in 
each grid evolves from the begin to the end of the simulation. Their kinetic energies are translated into 
the random thermal energetic particles by the "heating" processes in the downstream region resulting 
in a distinct enhancement temperature profiles in the shock front with the time. The profile of the ther- 
mal temperature shows the upstream averaged temperature of Tq = 2.5 x 10 5 A' and the downstream 
averaged temperature of Td = 9.0 x lO 6 ^. This means that the CME-driven shock can "heat" the 
background plasma efficiently and provide the first-order Fermi acceleration mechanism by crossing 
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the shock front for accelerating the particles, which are injected from the "heated" downstream region 
into the precursor region. 

Figure|2]shows a group of profiles of the physical parameters in the simulation. From top to bottom, 
the upper panel shows the proton density profile vs its position. The proton density is presented by 
the scaled value. The enhanced density flux apparently appears in the position of the shock front. The 
intensity of the density in the downstream is about five times larger than those in the unshocked upstream 
bulk flow. By comparing the proton density profiles between the the Figure [2] and the Figure [TJ the 
simulated bulk flow has a lightly higher proton density intensity in the downstream bulk flow than that 
in the observed downstream bulk flow. The middle panel in the Figure|2]denotes that the thermal velocity 
profile evolutes with the time. The profile at time of T~600 (it is zero before the simulation time T<600, 
take account of the injection from the PIB) has an initial Maxwellian thermal velocity of vq = 46km / s 
until it is shocked. After the profile is shocked as shown in middle panel of the Figure [2] it reaches an 
average thermal velocity with a value around < Vd >= 300fcm/s till the end of the simulation. Also 
the thermal velocity profile shows a slightly larger enhancement than that from the observation by Wind 
spacecraft. Although the simulated proton density and the solar wind thermal velocity are slightly larger 
than those from in-situ observations, we suggest that this is caused by the insufficiency of particles in 
the simulation. We have demonstrated that it is the case using a series of simulations with different 
initial number of particles per cell. The lower panel indicates the profiles of the bulk flow speed vs its 
position at the end of the simulation. The profile shows an upstream bulk flow speed U u = —QOOkm/s 
and the downstream bulk flow speed Ud = — 1042fcm/ s which are followed by the observations. The 
complex shock front fine structure will be showed in Figure [8] at the end of the simulation. The final 
evolutional positions of the FEB and the shock front are Xf = 31.5 and X s h = 121.5 in the x-axis, 
respectively. The distance between these two locations is just the size of the precursor region where the 
particle acceleration processes occur. It is just this region slowed the incoming upstream bulk flow speed 
U u down to the downstream bulk flow speed Ud- The bulk flow speed in precursor region is between the 
two bulk flow speeds (i.e. U u > U p > Ud)- From our simulation, we can see that the particle acceleration 
process and the "back pressure" due to the energetic particles occurred mostly in precursor region which 
results in a non-linear shock structure that is characterized by a bulk flow speed gradient. According to 
the evolutional shock front position X s h with the time, we can calculate the shock evolutional velocity 
V s h as follows. 

V\X max X s fo\ 
sh = = , (-5) 

-L max 

where, the X max is the total length of the simulation box, and the T max is the total simulation time. 
Then, we are able to calculate the total shock compression ratio in the shock frame as follows. 

? tot = , (4) 

\ V sh\ 

where AU is the relative bulk flow speed between the upstream and downstream, V s h is the shock 
velocity. 

Figure [8] shows the shock fine structure with the bulk flow speed near the shock front at the end 
of the simulation. V su b = 0.2103 shows the bulk flow speed of the subshock, Vd — shows the 
bulk flow speed of the downstream region, V s h = —0.0744 represents the value of the opposite shock 
evolutional velocity, and Uq = 0.3275 represents the incoming bulk flow speed with a related bulk flow 
speed difference of AU. All zones of the precursor, subshock and downstream are divided by a vertical 
dashed line and a solid line in the plot. These three zones constitute the total shock fine structure in 
the simulated shock region. The smooth precursor has a long scale in the range from the subshock's 
position X su i, to the FEB position Xf, which is invisible, beyond the left boundary of the plot. This 
zone is called the diffusive zone where the bulk flow speed will be slowed by the "back pressure" of 
the accelerated particles. The subshock region with a narrow scale of a three-grid-length has a deep 
drop of the bulk flow speed, in which the bulk flow speed vary from the subshock velocity V su b to the 
downstream velocity Vd- The scale of the three-grid-length is almost identical to the averaged thermal 
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Fig. 7 The mesh plot represents the evolutional bulk flow temperature profiles in their posi- 
tions with the time. The lower temperature represents the upstream bulk flow. The higher tem- 
perature represents the downstream bulk flow. The apparently distinguished boundary traces 
the shock front positions with time. 

mean free path over the downstream region. The subshock velocity V su b is decided by the horizontal 
dot-dashed line with a value of V su b = 0.2103. The downstream velocity Vd is marked with a horizontal 
dashed line at the end of the simulation, which should be with an averaged value of < Vd >= over 
the entire simulation time in the box frame. The nagative shock evolutional velocity marked with a 
horizontal solid line shows an value of V s h = —0.0744. We can calculate the subshock's compression 
ratio according to Rankie-Hongniout relationships in the shock frame as follows. 

Vsub + Vsh 

r sub — tt T/ 7 (•>) 

< Vd > +Vsh 

where, we take the averaged value of the downstream velocity < Vd > equal to zero. 
5 SUMMARY AND CONCLUSIONS 

In summary, we performed a dynamical Monte Carlo simulations on the 14-Dec-2006 CME-driven 
shock using an anisotropic scattering law. The specific temperature profile, shock fine structures, par- 
ticle injection function, as functions of time, are presented. We examined the correlation between the 
energy injection and the shock energy translation processes of the interplanetary CME-driven shock. 
Simultaneously, we find the simulated CME-driven shock energy spectrum provides a good fit to the 
observations from the multiple spacecraft. 

In conclusion, the dynamical Monte Carlo simulation of the 14-Dec-2006 CME-driven shock 
demonstrates that the energy spectrum is affected by the specific non-linear factor of the DSA. This 
paper focus on the energy injection, which is one of important nonlinear effects of the DSA. By cal- 
culating the energy injection rate of the CME-driven shock, we can understand how the CME-driven 
shock distributes its shock energy to accelerate the energetic particles by first-order Fermi acceleration 
mechanism as well as how it heats the solar wind background bulk flow at a certain efficiency. We give 
an energy injection rate of Ri n j = 15.32% in the 14-Dec-2006 CME-driven shock. We guess that this 
predicted injection rate could satisfy the required energies of the observed SEP events, which should be 
released from the CME-driven shock. 
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Fig. 8 The bulk flow speed fine structure of the simulated shock at the end of the simula- 
tion. The vertical dash line and vertical solid line split the entire region into three sections: 
precursor, subshock and downstream region. 
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